home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_300 / 334_03 / standard.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-02-04  |  15.5 KB  |  780 lines
  1. /* GNUPLOT - standard.c */
  2. /*
  3.  * Copyright (C) 1986, 1987, 1990   Thomas Williams, Colin Kelley
  4.  *
  5.  * Permission to use, copy, and distribute this software and its
  6.  * documentation for any purpose with or without fee is hereby granted, 
  7.  * provided that the above copyright notice appear in all copies and 
  8.  * that both that copyright notice and this permission notice appear 
  9.  * in supporting documentation.
  10.  *
  11.  * Permission to modify the software is granted, but not the right to
  12.  * distribute the modified code.  Modifications are to be distributed 
  13.  * as patches to released version.
  14.  *  
  15.  * This software  is provided "as is" without express or implied warranty.
  16.  * 
  17.  *
  18.  * AUTHORS
  19.  * 
  20.  *   Original Software:
  21.  *     Thomas Williams,  Colin Kelley.
  22.  * 
  23.  *   Gnuplot 2.0 additions:
  24.  *       Russell Lang, Dave Kotz, John Campbell.
  25.  * 
  26.  * send your comments or suggestions to (pixar!info-gnuplot@sun.com).
  27.  * 
  28.  */
  29.  
  30. #include <math.h>
  31. #include <stdio.h>
  32. #include "plot.h"
  33.  
  34. #ifdef vms
  35. #include <errno.h>
  36. #else
  37. extern int errno;
  38. #endif /* vms */
  39.  
  40.  
  41. extern struct value stack[STACK_DEPTH];
  42. extern int s_p;
  43. extern double zero;
  44.  
  45. struct value *pop(), *complex(), *integer();
  46.  
  47. double magnitude(), angle(), real(), imag();
  48.  
  49. /* The bessel function approximations here are from
  50.  * "Computer Approximations"
  51.  * by Hart, Cheney et al.
  52.  * John Wiley & Sons, 1968
  53.  */
  54.  
  55. /* There appears to be a mistake in Hart, Cheney et al. on page 149.
  56.  * Where it list Qn(x)/x ~ P(z*z)/Q(z*z), z = 8/x, it should read
  57.  *               Qn(x)/z ~ P(z*z)/Q(z*z), z = 8/x
  58.  * In the functions below, Qn(x) is implementated using the later
  59.  * equation.
  60.  * These bessel functions are accurate to about 1e-13
  61.  */
  62.  
  63. #define PI_ON_FOUR       0.78539816339744830961566084581987572
  64. #define PI_ON_TWO        1.57079632679489661923131269163975144
  65. #define THREE_PI_ON_FOUR 2.35619449019234492884698253745962716
  66. #define TWO_ON_PI        0.63661977236758134307553505349005744
  67.  
  68. static double dzero = 0.0;
  69.  
  70. /* jzero for x in [0,8]
  71.  * Index 5849, 19.22 digits precision
  72.  */
  73. static double pjzero[] = {
  74.      0.4933787251794133561816813446e+21,
  75.     -0.11791576291076105360384408e+21,
  76.      0.6382059341072356562289432465e+19,
  77.     -0.1367620353088171386865416609e+18,
  78.      0.1434354939140346111664316553e+16,
  79.     -0.8085222034853793871199468171e+13,
  80.      0.2507158285536881945555156435e+11,
  81.     -0.4050412371833132706360663322e+8,
  82.      0.2685786856980014981415848441e+5
  83. };
  84.  
  85. static double qjzero[] = {
  86.     0.4933787251794133562113278438e+21,
  87.     0.5428918384092285160200195092e+19,
  88.     0.3024635616709462698627330784e+17,
  89.     0.1127756739679798507056031594e+15,
  90.     0.3123043114941213172572469442e+12,
  91.     0.669998767298223967181402866e+9,
  92.     0.1114636098462985378182402543e+7,
  93.     0.1363063652328970604442810507e+4,
  94.     0.1e+1
  95. };
  96.  
  97. /* pzero for x in [8,inf]
  98.  * Index 6548, 18.16 digits precision
  99.  */
  100. static double ppzero[] = {
  101.     0.2277909019730468430227002627e+5,
  102.     0.4134538663958076579678016384e+5,
  103.     0.2117052338086494432193395727e+5,
  104.     0.348064864432492703474453111e+4,
  105.     0.15376201909008354295771715e+3,
  106.     0.889615484242104552360748e+0
  107. };
  108.  
  109. static double qpzero[] = {
  110.     0.2277909019730468431768423768e+5,
  111.     0.4137041249551041663989198384e+5,
  112.     0.2121535056188011573042256764e+5,
  113.     0.350287351382356082073561423e+4,
  114.     0.15711159858080893649068482e+3,
  115.     0.1e+1
  116. };
  117.  
  118. /* qzero for x in [8,inf]
  119.  * Index 6948, 18.33 digits precision
  120.  */
  121. static double pqzero[] = {
  122.     -0.8922660020080009409846916e+2,
  123.     -0.18591953644342993800252169e+3,
  124.     -0.11183429920482737611262123e+3,
  125.     -0.2230026166621419847169915e+2,
  126.     -0.124410267458356384591379e+1,
  127.     -0.8803330304868075181663e-2,
  128. };
  129.  
  130. static double qqzero[] = {
  131.     0.571050241285120619052476459e+4,
  132.     0.1195113154343461364695265329e+5,
  133.     0.726427801692110188369134506e+4,
  134.     0.148872312322837565816134698e+4,
  135.     0.9059376959499312585881878e+2,
  136.     0.1e+1
  137. };
  138.  
  139.  
  140. /* yzero for x in [0,8]
  141.  * Index 6245, 18.78 digits precision
  142.  */
  143. static double pyzero[] = {
  144.     -0.2750286678629109583701933175e+20,
  145.      0.6587473275719554925999402049e+20,
  146.     -0.5247065581112764941297350814e+19,
  147.      0.1375624316399344078571335453e+18,
  148.     -0.1648605817185729473122082537e+16,
  149.      0.1025520859686394284509167421e+14,
  150.     -0.3436371222979040378171030138e+11,
  151.      0.5915213465686889654273830069e+8,
  152.     -0.4137035497933148554125235152e+5
  153. };
  154.  
  155. static double qyzero[] = {
  156.     0.3726458838986165881989980739e+21,
  157.     0.4192417043410839973904769661e+19,
  158.     0.2392883043499781857439356652e+17,
  159.     0.9162038034075185262489147968e+14,
  160.     0.2613065755041081249568482092e+12,
  161.     0.5795122640700729537380087915e+9,
  162.     0.1001702641288906265666651753e+7,
  163.     0.1282452772478993804176329391e+4,
  164.     0.1e+1
  165. };
  166.  
  167.  
  168. /* jone for x in [0,8]
  169.  * Index 6050, 20.98 digits precision
  170.  */
  171. static double pjone[] = {
  172.      0.581199354001606143928050809e+21,
  173.     -0.6672106568924916298020941484e+20,
  174.      0.2316433580634002297931815435e+19,
  175.     -0.3588817569910106050743641413e+17,
  176.      0.2908795263834775409737601689e+15,
  177.     -0.1322983480332126453125473247e+13,
  178.      0.3413234182301700539091292655e+10,
  179.     -0.4695753530642995859767162166e+7,
  180.      0.270112271089232341485679099e+4
  181. };
  182.  
  183. static double qjone[] = {
  184.     0.11623987080032122878585294e+22,
  185.     0.1185770712190320999837113348e+20,
  186.     0.6092061398917521746105196863e+17,
  187.     0.2081661221307607351240184229e+15,
  188.     0.5243710262167649715406728642e+12,
  189.     0.1013863514358673989967045588e+10,
  190.     0.1501793594998585505921097578e+7,
  191.     0.1606931573481487801970916749e+4,
  192.     0.1e+1
  193. };
  194.  
  195.  
  196. /* pone for x in [8,inf]
  197.  * Index 6749, 18.11 digits precision
  198.  */
  199. static double ppone[] = {
  200.     0.352246649133679798341724373e+5,
  201.     0.62758845247161281269005675e+5,
  202.     0.313539631109159574238669888e+5,
  203.     0.49854832060594338434500455e+4,
  204.     0.2111529182853962382105718e+3,
  205.     0.12571716929145341558495e+1
  206. };
  207.  
  208. static double qpone[] = {
  209.     0.352246649133679798068390431e+5,
  210.     0.626943469593560511888833731e+5,
  211.     0.312404063819041039923015703e+5,
  212.     0.4930396490181088979386097e+4,
  213.     0.2030775189134759322293574e+3,
  214.     0.1e+1
  215. };
  216.  
  217. /* qone for x in [8,inf]
  218.  * Index 7149, 18.28 digits precision
  219.  */
  220. static double pqone[] = {
  221.     0.3511751914303552822533318e+3,
  222.     0.7210391804904475039280863e+3,
  223.     0.4259873011654442389886993e+3,
  224.     0.831898957673850827325226e+2,
  225.     0.45681716295512267064405e+1,
  226.     0.3532840052740123642735e-1
  227. };
  228.  
  229. static double qqone[] = {
  230.     0.74917374171809127714519505e+4,
  231.     0.154141773392650970499848051e+5,
  232.     0.91522317015169922705904727e+4,
  233.     0.18111867005523513506724158e+4,
  234.     0.1038187585462133728776636e+3,
  235.     0.1e+1
  236. };
  237.  
  238.  
  239. /* yone for x in [0,8]
  240.  * Index 6444, 18.24 digits precision
  241.  */
  242. static double pyone[] = {
  243.     -0.2923821961532962543101048748e+20,
  244.      0.7748520682186839645088094202e+19,
  245.     -0.3441048063084114446185461344e+18,
  246.      0.5915160760490070618496315281e+16,
  247.     -0.4863316942567175074828129117e+14,
  248.      0.2049696673745662182619800495e+12,
  249.     -0.4289471968855248801821819588e+9,
  250.      0.3556924009830526056691325215e+6
  251. };
  252.  
  253. static double qyone[] = {
  254.     0.1491311511302920350174081355e+21,
  255.     0.1818662841706134986885065935e+19,
  256.     0.113163938269888452690508283e+17,
  257.     0.4755173588888137713092774006e+14,
  258.     0.1500221699156708987166369115e+12,
  259.     0.3716660798621930285596927703e+9,
  260.     0.726914730719888456980191315e+6,
  261.     0.10726961437789255233221267e+4,
  262.     0.1e+1
  263. };
  264.  
  265.  
  266. f_real()
  267. {
  268. struct value a;
  269.     push( complex(&a,real(pop(&a)), 0.0) );
  270. }
  271.  
  272. f_imag()
  273. {
  274. struct value a;
  275.     push( complex(&a,imag(pop(&a)), 0.0) );
  276. }
  277.  
  278. f_arg()
  279. {
  280. struct value a;
  281.     push( complex(&a,angle(pop(&a)), 0.0) );
  282. }
  283.  
  284. f_conjg()
  285. {
  286. struct value a;
  287.     (void) pop(&a);
  288.     push( complex(&a,real(&a),-imag(&a) ));
  289. }
  290.  
  291. f_sin()
  292. {
  293. struct value a;
  294.     (void) pop(&a);
  295.     push( complex(&a,sin(real(&a))*cosh(imag(&a)), cos(real(&a))*sinh(imag(&a))) );
  296. }
  297.  
  298. f_cos()
  299. {
  300. struct value a;
  301.     (void) pop(&a);
  302.     push( complex(&a,cos(real(&a))*cosh(imag(&a)), -sin(real(&a))*sinh(imag(&a))));
  303. }
  304.  
  305. f_tan()
  306. {
  307. struct value a;
  308. register double den;
  309.     (void) pop(&a);
  310.     if (imag(&a) == 0.0)
  311.         push( complex(&a,tan(real(&a)),0.0) );
  312.     else {
  313.         den = cos(2*real(&a))+cosh(2*imag(&a));
  314.         if (den == 0.0) {
  315.             undefined = TRUE;
  316.             push( &a );
  317.         }
  318.         else
  319.             push( complex(&a,sin(2*real(&a))/den, sinh(2*imag(&a))/den) );
  320.     }
  321. }
  322.  
  323. f_asin()
  324. {
  325. struct value a;
  326. register double alpha, beta, x, y;
  327.     (void) pop(&a);
  328.     x = real(&a); y = imag(&a);
  329.     if (y == 0.0) {
  330.         if (fabs(x) > 1.0) {
  331.             undefined = TRUE;
  332.             push(complex(&a,0.0, 0.0));
  333.         } else
  334.             push( complex(&a,asin(x),0.0) );
  335.     } else {
  336.         beta  = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
  337.         alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
  338.         push( complex(&a,asin(beta), log(alpha + sqrt(alpha*alpha-1))) );
  339.     }
  340. }
  341.  
  342. f_acos()
  343. {
  344. struct value a;
  345. register double alpha, beta, x, y;
  346.     (void) pop(&a);
  347.     x = real(&a); y = imag(&a);
  348.     if (y == 0.0) {
  349.         if (fabs(x) > 1.0) {
  350.             undefined = TRUE;
  351.             push(complex(&a,0.0, 0.0));
  352.         } else
  353.             push( complex(&a,acos(x),0.0) );
  354.     } else {
  355.         alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
  356.         beta  = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
  357.         push( complex(&a,acos(beta), log(alpha + sqrt(alpha*alpha-1))) );
  358.     }
  359. }
  360.  
  361. f_atan()
  362. {
  363. struct value a;
  364. register double x, y;
  365.     (void) pop(&a);
  366.     x = real(&a); y = imag(&a);
  367.     if (y == 0.0)
  368.         push( complex(&a,atan(x), 0.0) );
  369.     else if (x == 0.0 && fabs(y) == 1.0) {
  370.         undefined = TRUE;
  371.         push(complex(&a,0.0, 0.0));
  372.     } else
  373.         push( complex(&a,atan(2*x/(1-x*x-y*y)),
  374.                 log((x*x+(y+1)*(y+1))/(x*x+(y-1)*(y-1)))/4) );
  375. }
  376.  
  377. f_sinh()
  378. {
  379. struct value a;
  380.     (void) pop(&a);
  381.     push( complex(&a,sinh(real(&a))*cos(imag(&a)), cosh(real(&a))*sin(imag(&a))) );
  382. }
  383.  
  384. f_cosh()
  385. {
  386. struct value a;
  387.     (void) pop(&a);
  388.     push( complex(&a,cosh(real(&a))*cos(imag(&a)), sinh(real(&a))*sin(imag(&a))) );
  389. }
  390.  
  391. f_tanh()
  392. {
  393. struct value a;
  394. register double den;
  395.     (void) pop(&a);
  396.     den = cosh(2*real(&a)) + cos(2*imag(&a));
  397.     push( complex(&a,sinh(2*real(&a))/den, sin(2*imag(&a))/den) );
  398. }
  399.  
  400. f_int()
  401. {
  402. struct value a;
  403.     push( integer(&a,(int)real(pop(&a))) );
  404. }
  405.  
  406.  
  407. f_abs()
  408. {
  409. struct value a;
  410.     (void) pop(&a);
  411.     switch (a.type) {
  412.         case INT:
  413.             push( integer(&a,abs(a.v.int_val)) );            
  414.             break;
  415.         case CMPLX:
  416.             push( complex(&a,magnitude(&a), 0.0) );
  417.     }
  418. }
  419.  
  420. f_sgn()
  421. {
  422. struct value a;
  423.     (void) pop(&a);
  424.     switch(a.type) {
  425.         case INT:
  426.             push( integer(&a,(a.v.int_val > 0) ? 1 : 
  427.                     (a.v.int_val < 0) ? -1 : 0) );
  428.             break;
  429.         case CMPLX:
  430.             push( integer(&a,(a.v.cmplx_val.real > 0.0) ? 1 : 
  431.                     (a.v.cmplx_val.real < 0.0) ? -1 : 0) );
  432.             break;
  433.     }
  434. }
  435.  
  436.  
  437. f_sqrt()
  438. {
  439. struct value a;
  440. register double mag, ang;
  441.     (void) pop(&a);
  442.     mag = sqrt(magnitude(&a));
  443.     if (imag(&a) == 0.0 && real(&a) < 0.0)
  444.         push( complex(&a,0.0,mag) );
  445.     else
  446.     {
  447.         if ( (ang = angle(&a)) < 0.0)
  448.             ang += 2*Pi;
  449.         ang /= 2;
  450.         push( complex(&a,mag*cos(ang), mag*sin(ang)) );
  451.     }
  452. }
  453.  
  454.  
  455. f_exp()
  456. {
  457. struct value a;
  458. register double mag, ang;
  459.     (void) pop(&a);
  460.     mag = exp(real(&a));
  461.     ang = imag(&a);
  462.     push( complex(&a,mag*cos(ang), mag*sin(ang)) );
  463. }
  464.  
  465.  
  466. f_log10()
  467. {
  468. struct value a;
  469. register double l10;;
  470.     (void) pop(&a);
  471.     l10 = log(10.0);    /***** replace with a constant! ******/
  472.     push( complex(&a,log(magnitude(&a))/l10, angle(&a)/l10) );
  473. }
  474.  
  475.  
  476. f_log()
  477. {
  478. struct value a;
  479.     (void) pop(&a);
  480.     push( complex(&a,log(magnitude(&a)), angle(&a)) );
  481. }
  482.  
  483.  
  484. f_floor()
  485. {
  486. struct value a;
  487.  
  488.     (void) pop(&a);
  489.     switch (a.type) {
  490.         case INT:
  491.             push( integer(&a,(int)floor((double)a.v.int_val)));            
  492.             break;
  493.         case CMPLX:
  494.             push( integer(&a,(int)floor(a.v.cmplx_val.real)));
  495.     }
  496. }
  497.  
  498.  
  499. f_ceil()
  500. {
  501. struct value a;
  502.  
  503.     (void) pop(&a);
  504.     switch (a.type) {
  505.         case INT:
  506.             push( integer(&a,(int)ceil((double)a.v.int_val)));            
  507.             break;
  508.         case CMPLX:
  509.             push( integer(&a,(int)ceil(a.v.cmplx_val.real)));
  510.     }
  511. }
  512.  
  513. #ifdef GAMMA
  514.  
  515. f_gamma()
  516. {
  517. extern int signgam;
  518. register double y;
  519. struct value a;
  520.  
  521.     y = GAMMA(real(pop(&a)));
  522.     if (y > 88.0) {
  523.         undefined = TRUE;
  524.         push( integer(&a,0) );
  525.     }
  526.     else
  527.         push( complex(&a,signgam * exp(y),0.0) );
  528. }
  529.  
  530. #endif /* GAMMA */
  531.  
  532.  
  533. /* bessel function approximations */
  534. double jzero(x)
  535. double x;
  536. {
  537. double p, q, x2;
  538. int n;
  539.  
  540.     x2 = x * x;
  541.     p = pjzero[8];
  542.     q = qjzero[8];
  543.     for (n=7; n>=0; n--) {
  544.         p = p*x2 + pjzero[n];
  545.         q = q*x2 + qjzero[n];
  546.     }
  547.     return(p/q);
  548. }
  549.  
  550. double pzero(x)
  551. double x;
  552. {
  553. double p, q, z, z2;
  554. int n;
  555.  
  556.     z = 8.0 / x;
  557.     z2 = z * z;
  558.     p = ppzero[5];
  559.     q = qpzero[5];
  560.     for (n=4; n>=0; n--) {
  561.         p = p*z2 + ppzero[n];
  562.         q = q*z2 + qpzero[n];
  563.     }
  564.     return(p/q);
  565. }
  566.  
  567. double qzero(x)
  568. double x;
  569. {
  570. double p, q, z, z2;
  571. int n;
  572.  
  573.     z = 8.0 / x;
  574.     z2 = z * z;
  575.     p = pqzero[5];
  576.     q = qqzero[5];
  577.     for (n=4; n>=0; n--) {
  578.         p = p*z2 + pqzero[n];
  579.         q = q*z2 + qqzero[n];
  580.     }
  581.     return(p/q);
  582. }
  583.  
  584. double yzero(x)
  585. double x;
  586. {
  587. double p, q, x2;
  588. int n;
  589.  
  590.     x2 = x * x;
  591.     p = pyzero[8];
  592.     q = qyzero[8];
  593.     for (n=7; n>=0; n--) {
  594.         p = p*x2 + pyzero[n];
  595.         q = q*x2 + qyzero[n];
  596.     }
  597.     return(p/q);
  598. }
  599.  
  600. double rj0(x)
  601. double x;
  602. {
  603.     if ( x <= 0.0 )
  604.         x = -x;
  605.     if ( x < 8.0 )
  606.         return(jzero(x));
  607.     else
  608.         return( sqrt(TWO_ON_PI/x) *
  609.             (pzero(x)*cos(x-PI_ON_FOUR) - 8.0/x*qzero(x)*sin(x-PI_ON_FOUR)) );
  610.  
  611. }
  612.  
  613. double ry0(x)
  614. double x;
  615. {
  616.     if ( x < 0.0 )
  617.         return(dzero/dzero); /* error */
  618.     if ( x < 8.0 )
  619.         return( yzero(x) + TWO_ON_PI*rj0(x)*log(x) );
  620.     else
  621.         return( sqrt(TWO_ON_PI/x) *
  622.             (pzero(x)*sin(x-PI_ON_FOUR) + 
  623.             (8.0/x)*qzero(x)*cos(x-PI_ON_FOUR)) );
  624.  
  625. }
  626.  
  627.  
  628. double jone(x)
  629. double x;
  630. {
  631. double p, q, x2;
  632. int n;
  633.  
  634.     x2 = x * x;
  635.     p = pjone[8];
  636.     q = qjone[8];
  637.     for (n=7; n>=0; n--) {
  638.         p = p*x2 + pjone[n];
  639.         q = q*x2 + qjone[n];
  640.     }
  641.     return(p/q);
  642. }
  643.  
  644. double pone(x)
  645. double x;
  646. {
  647. double p, q, z, z2;
  648. int n;
  649.  
  650.     z = 8.0 / x;
  651.     z2 = z * z;
  652.     p = ppone[5];
  653.     q = qpone[5];
  654.     for (n=4; n>=0; n--) {
  655.         p = p*z2 + ppone[n];
  656.         q = q*z2 + qpone[n];
  657.     }
  658.     return(p/q);
  659. }
  660.  
  661. double qone(x)
  662. double x;
  663. {
  664. double p, q, z, z2;
  665. int n;
  666.  
  667.     z = 8.0 / x;
  668.     z2 = z * z;
  669.     p = pqone[5];
  670.     q = qqone[5];
  671.     for (n=4; n>=0; n--) {
  672.         p = p*z2 + pqone[n];
  673.         q = q*z2 + qqone[n];
  674.     }
  675.     return(p/q);
  676. }
  677.  
  678. double yone(x)
  679. double x;
  680. {
  681. double p, q, x2;
  682. int n;
  683.  
  684.     x2 = x * x;
  685.     p = 0.0;
  686.     q = qyone[8];
  687.     for (n=7; n>=0; n--) {
  688.         p = p*x2 + pyone[n];
  689.         q = q*x2 + qyone[n];
  690.     }
  691.     return(p/q);
  692. }
  693.  
  694. double rj1(x)
  695. double x;
  696. {
  697. double v,w;
  698.     v = x;
  699.     if ( x < 0.0 )
  700.         x = -x;
  701.     if ( x < 8.0 )
  702.         return(v*jone(x));
  703.     else {
  704.         w = sqrt(TWO_ON_PI/x) *
  705.             (pone(x)*cos(x-THREE_PI_ON_FOUR) - 
  706.                8.0/x*qone(x)*sin(x-THREE_PI_ON_FOUR)) ;
  707.         if (v < 0.0)
  708.             w = -w;
  709.         return( w );
  710.     }
  711. }
  712.  
  713. double ry1(x)
  714. double x;
  715. {
  716.     if ( x <= 0.0 )
  717.         return(dzero/dzero); /* error */
  718.     if ( x < 8.0 )
  719.         return( x*yone(x) + TWO_ON_PI*(rj1(x)*log(x) - 1.0/x) );
  720.     else
  721.         return( sqrt(TWO_ON_PI/x) *
  722.             (pone(x)*sin(x-THREE_PI_ON_FOUR) + 
  723.             (8.0/x)*qone(x)*cos(x-THREE_PI_ON_FOUR)) );
  724. }
  725.  
  726.  
  727. f_besj0()    
  728. {
  729. struct value a;
  730. double x;
  731.     (void) pop(&a);
  732.     if (imag(&a) > zero)
  733.         int_error("can only do bessel functions of reals",NO_CARET);
  734.     push( complex(&a,rj0(real(&a)),0.0) );
  735. }
  736.  
  737.  
  738. f_besj1()    
  739. {
  740. struct value a;
  741. double x;
  742.     (void) pop(&a);
  743.     if (imag(&a) > zero)
  744.         int_error("can only do bessel functions of reals",NO_CARET);
  745.     push( complex(&a,rj1(real(&a)),0.0) );
  746. }
  747.  
  748.  
  749. f_besy0()    
  750. {
  751. struct value a;
  752. double x;
  753.     (void) pop(&a);
  754.     if (imag(&a) > zero)
  755.         int_error("can only do bessel functions of reals",NO_CARET);
  756.     if (real(&a) > 0.0)
  757.         push( complex(&a,ry0(real(&a)),0.0) );
  758.     else {
  759.         push( complex(&a,0.0,0.0) );
  760.         undefined = TRUE ;
  761.     }
  762. }
  763.  
  764.  
  765. f_besy1()    
  766. {
  767. struct value a;
  768. double x;
  769.     (void) pop(&a);
  770.     if (imag(&a) > zero)
  771.         int_error("can only do bessel functions of reals",NO_CARET);
  772.     if (real(&a) > 0.0)
  773.         push( complex(&a,ry1(real(&a)),0.0) );
  774.     else {
  775.         push( complex(&a,0.0,0.0) );
  776.         undefined = TRUE ;
  777.     }
  778. }
  779.  
  780.